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Notation 


a auxiliary variable 

b wing span 

B range, force of inertia 

c constant 

C solid surface 

D image determinant, "thickness" 

e Importance function 

e viscosity parameter, small sphere 

r\ image coordinate 

f function 

4» potential 

4 disturbance potential 

G form function 

r circulation 

H hyperbolic (auxiliary) point 

1 counting index 

K wake 

K adiabatic exponent 

L corridor 

M Mach number 

m mass, "middle" 

n edge perpendicular 

0 surface element 

o "top", "zero" and "dynamic air pressure" 

P pressure force, check point 

p pressure 

q total velocity 

r radius 

p density 

s path distance 

u,v,w velocity components 
X coordinate 

X auxiliary function 


K 

image coordinate 


y 

coordinate 


Y 

auxiliary function 


z 

coordinate 


Z 

auxiliary function 


c 

image coordinate 


« 

"sound” 


List of PlKures 


Pig. 


Page No 

1 

NACA 0012-wing, input 

24 

2-16 

Pressure distributions 

25-39 

17 

Coefficients 

40 

18-19 

Pressure distribution relief 

41-42 

20-21 

Coefficient distributions 

43-44 

22 

Center of pressure 

45 

23 

BlOX-wing, input 

46 

24-33 

Pressure distributions 

47-56 

34 

Coefficients 

57 

35-36 

Pressure distribution relief 

58-59 

37-39 

Coefficient distributions 

60-62 

40 

Center of pressure 

63 


Hi 


Table of Contents Page 

1. Introduction 1 

2. Fundamental Equations 1 

3i Standardization 2 

Distant Field Solution 3 

4.1 Prandtl Transformation 3 

4.2 Integral Potential Representation 5 

4.3 Displacement Term 7 

4.4 Vortex Term 8 

5. Continuity Equation Modification for 

Transonic Flow 9 

5.1 Directional Dependence in the Case of 

Supersonic Plow 9 

5.2 Artificial Viscosity 10 

5.3 Viscosity Parameter 11 

6. The Variation Principle 12 

6.1 The Euler Minimum Principle 12 

6.2 The Pressure Minimum Integral 13 

7. Numerical Evaluation 14 

7.1 The Eight-Node Element Cube l4 

7.2 The General Eight-Node Volume Element 15 

7.3 The Element Pressure Minimum Integral 17 

7.4 The Global Pressure Minimum Integral l8 

7.5 Relaxation Method l8 

8. Results 20 

9. Summary 22 


iv 


A FINITE VOLUME METHOD FOR CALCULATING TRANSONIC POTENTIAL FLOW 
AROUND WINGS FROM THE PRESSURE MINIMUM INTEGRAL 


Albrecht Eberle 

Messerschmltt-Bdlkow-Blohm GmbH, Ottorbrunn near Munich, 

West Germany 


1,. Introduction /I 

The development of a successful finite element method for 
calculating transonic flow around a profile [1] provided the 
precondition for progranunlng this method. 

In contrast to the usual difference methods, finite volume 
methods do not require a rectangular mesh network so that they are 
particularly suitable for treating complex aerodynamic configur- 
ations. Even if this paper only covers airfoil wing calculations 
In the case of transonic oncoming flow, naturally one has in the 
back of one’s mind the Intention of expanding the computer 
program in the above direction when the occasion arises. 

The introduction of an extremely simply formulated concept 
of plastic viscosity makes it possible to calculate shock- 
effected supercritical pressure distributions on any wings. 

If, moreover, the distant field solution for supersonic 
flow is used, this offers the possibility of calculating pressure 
distributions for oncoming flow mach numbers even greater than 
one. 

2. Fundamental Equations /2 

If we assume frictionless stationary flow of an ideal gas, 
this can be described by the continuity equation 

(^u)^ + ® 

» Numbers in the margin indicate pagination in the foreign text. 
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and the condition of Irrotatlonallty 


Uz 

Uw 



CVJ 

• Vx 

(3) 

“ Wy 

(^) 


With uniform oncoming flow the energy carried with the fluid 
element is constant: 


2tu Po 

*-i s ^ ^ JZ 


(5) 


with 


5 O ? 

+ w 


The change of state is determined by the following isentrope 




with the Mach number M defined as follows 


( 6 ) 


% 


_E. = ^ 


8 

then by inserting this in Eq. (5) the squared speed becomes 

XPo 


>2 « 


8o 


(7) 


3. Standardization 


The Irrelevant constants in Eqs. (5), (6) and 
eliminated by introducing a dimensionless velocity 


(7) can be 

q: 


q 


2 




/3 


2 


Then Eq. (5) becomes 


9oP -2 
Po8 ^ 


1 


and Eq. (7) becomes 


-2 

q 



M 


2 


( 8 ) 


(9) 


Eq. (6) In (8) gives us: 

__2 

^ (1 - q )x-l (10) 

Prom Eq. (9) we read the standardized sound velocity as 



If the undeleted variables are now understood to be the 
appropriately standardized variables, equations (1), (2), 

(3) and (^) remain unchanged. 

Distant Field Solution 
^.1 Prandtl Transformation 

The flow behavior far downstream from the wing can be 

determined by the approximate solution of Eq. (1) for 
2 2 2 

X + y + z In this regard, Eq. (1) is differentiated out 

and in so doing the transverse velocity components with respect 
to u are Ignored. 


uOv + Vy + Wg) “ 0 

^ ^ ( 11 ) 

the term is determined by Eq. (10) with back differentiation: 


2g 

(x-t) (l-q2) 




■’Sq2 ■ 


The squared velocity is replaced by Eq. (9); 

Ux 

If we make the following approximate assumption that M “ then 
Eq. (11) becomes: 


(1 - m2) Ux * Vy * Wj . 0 

Ac a secondary condition, Eqs. (2), (3) and (M) have to be 
fulfilled. This a accomplished by Introducing a potential 
according to the following specification: 


= <t>x- V • <t>y, W » (j)j^ 


(13) 


as a result of which Eq. (12) assumes the following form; 

(1 - m£) ())xx ^ 4>yY ^ ° 

If we transform the y and z axes with 


1) = "I 1 - Y, E = 1 - z 

then the Laplace equation becomes 


(IM 




4^xx ♦ 


0 


( 15 ) 


A. 2 Integral Potential Representation 
If we insert the equation 


Z5 





(16) 


for a turbulent parallel main flow into Eq. (15), then obviously 
the following equation must be solved: 


“Pxx ‘‘Ptiti - 0 

To this end, Eq. (17) is multiplied by an importance function 
e and partially integrated over the region of flow, taking Into 
consideration the boundary shown in the sketch below: 



e (^>xx 

-KLE ^ 

- ^ e<pn dO dxdT|d^ 

CK^LE (B)-KLE 


ni 


0 


( 18 ) 
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If we replace e with we similarly obtain; 


®xx dxdTjd^ » 

(B)-KL£ 

* " JJJ ^®xVx ■*■ ®7} '^7]+ dxdT)d^ 

CKe« L£ (B)-KLE 

The volume Integral on the right side Is replaced by Eq. (l8): 

§ ('Psn - eiPn> <10 " JJJiplexx * em * err) <Jxdi)dr 

CK.LE (B)-KUe 

According to the definition and taking into consideration Eq, (16) 
the normal derivative and the disturbance ootential ^ disappear 
at infinity. 

On the corridor L the derivatives for outgoing flight and 
return flight cancel each other out so that their contribution 
disappears . 

In the wake K, the normal derivatives of the disturbance 
potential Indeed cancel each other out, since there is a steady 
velocity field, but a potential difference between the top and 
bottom side should be allowed in order to simulate the departing 
vortex layer. This thus leaves us with the following equation 
from the Integral equation; 

^ipe^dO = ^eipndO + § (eq)^, 

EEC 

+ JJJ^^exx * eyy + 

(B)-£K 


-<pep) dO -J|«4>end0 
K 

dxd*qd^ 


The left Integral can be evaluated Immediately if e represents 

the surface of a small sphere and e Is exclusively a function of 

the radius. Then ♦ tends towards the value and 6/6r. becomes 

P 

6/6r, with the Integrand becoming a constant; 

^>€^<30 = ^dO = 4Tlr^ 

G 


§ 


If we now call for 


Oj- 4 TC = 1 


It follows that 


1 

® = “ TtTF 


If we Insert this lesult Into the integral equation and let 
E shrink to 0, we finally obtain: 


=■ '*0 '’ o ] 

C + K C 


(19) 


^.3 Displacement Term /7 

We are looking for an approximate solution of the second 
Integral of Eq. (19) for a wing of moderate thickness. Then the 
following simplifications are valid; 

dO » dxdT] » dx - m 2 dy 

^ , 3 

If 1 - d z (from (1^)) 

C 

For small excess velocities the linearized limit condition 


7 

I 


can be used, with which the partial integral can be Integrated 
as follows; 






rear 

front 




The first term does not apply in the case of closed profile 
sections. Also, far downstream from the wing the effect of 
singularities appears punctiform so that ultimately we arrive 
at the following final result: 




^op r ~ 

^ [(yp-y^)^ + (zp-z„)2] I 3/2 




zdx • dy 
section 


( 20 ) 


The integrations are performed numerically in section. 
Vortex Term 

The first Integral from Eq. (19) can be evaluated as 
follows using appropriate linearizations: 



C + K 


(continued on next page) 


b 00 

^ P f r A <4? dxdy => 

4TI J J r3 

-bo 


A^)d (xp-x) dy 


_ _iL f r 

-i i {(xp-x) 2 *(l-Mj,)[(zp-r) 2 .(yp-y) 2 ] j 3/2 

b oo 

2 p r r A^>[a 2 + (xp-x)2_(xp-x)2 ] d(xp-x)dy 
i i [a2 . (xp-x>n 3/2 


zp p A^ (Xp-x) 
^J72"Va"2 V (xp-x)J' 


CO 

dy 

m 


In these equations, because of the constancy of the velocities 
on the top and underside of the wake, the potential Jump was 
assumed constant. We thus obtain: 



Zp 

4Tl 


J 



Avp 

+ (yp-y)^ 



Xp-x 


m 


|(xp-Xn^)2+(l_M^) [ (yp-y) 2+ (zp_ 2 ^) 2 ]' j 


dy 


( 21 ) 


where the m line can be assumed to be at about 25 % of the local 
chord length. 


5» Continuity Equation Modification for Transonic Flow /9 

5.1 Directional Dependence in the Case of Supersonic Flow 

Without a provision to take into consideration the change 
in the type of flow at transonic velocities, a numerical algorithm 


9 


for solving continuity equation (1) would fall. Moreover, shock 
waves could not develop automatically In the course of the 
calculation. If we now want to determine the velocity potential 
at a special point P of the supersonic region, it must be borne 
in mind that the chebk point is influenced exclusively by physical 
data from the accompanying forward Mach cone. 

It therefore suggests itself to coordinate the values pu, pv 
and pw from Eq. (1) at a point H, which lies a small distance 
upstream from P, to the check point itself. This ensures that 
signals picked up downstream cannot reach the check point . In 
the case of subsonic flow, on the other hand, there is not 
directional dependency. Here the physical variables at the check 
point itself are used in the calculation. Accordingly, in the 
case of transonic flow, a numerical case distinction must be 
performed which depends on whether local supersonic or subsonic 
velocity is present. 

5.2 Artificial Viscosity 

The suggested assignment of physical quantities to check 
point P in the case of supersonic flow can be formulated 
mathematically by the linear upstream development of the variables 
pu, pv and pw from Eq. (1). This should first of all be done for 
pu: 

(9u)h = (§u)p + (su)ps - Sp) 


where S is at least approximately the run length along the flow 
line through P. The derivation (pu)g is formed as follows: 

^(b5) ^ b ^ ^ b ^ S(nS) 

Since H is only a small distance from P, the flow line can be 



approximated as a straight line so that (u/q)g disappears. 

(^u)^ = ^ (Sq)s ^ 2 

Since the local density must be computed In the computer pro- 
gram anyway. It Is advisable to change over from S to p as the 
Independent variable : 


(^q)s As = Ag 

Differentiation with the help of Eqa. (10) and (9) gives us 

(§q)sAS = q (1 - ^)Ag . 

For distinguishing the type of flow we can here comfortably 
Introduce the switching function "max": 

(gu)n = u [g + max (0, 1 - <§h 


For the terms pv, pw It Is only necessary to replace u with v. 

A transonic computing process thus works If the following 
expression Is merely placed In the check point for the density: 

j — ► J + max (0, 1 - ) A ^ 


5.3 Viscosity Parameter /H 

Since the vector length between auxiliary point H and check 
point P can be chosen more or less arbitrarily, the artificial 
viscosity Is generalized by a parameter e which Is constant In the 
entire flow field; 

^ + e max (0, 1 - A^ 
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Thus, on the one hand, the accuracy of the method can be 
Increased (small s), and on the other hand iterative convergence 
from case to case can be guaranteed (large e). 

6, The Variation Principle /12 

6.1 The Euler Minimum Principle 

Here we will try to put the continuity equation (1) in a 
form appropriate for the computer which does not allow the 
introduction of a nonperpendicular calculating mesh network. 

In this regard, in reference [1] a variation principle was 
derived from the weighted remainder and in reference [2] from 
the least square. An especially illuminating procedure is based 
on the Euler principle according to which the equilibrium of 
forces is formed on the fluid particle. 



This should be done here first of all for the x-dlrectlon: 


pdydz - mu - (p + Pxdx) dydz = 0 
Px • = 0 

In the case of stationary flow, the acceleration u can be 
written as follows independent of time: 

u » u^x + Uyy + u^z « uux vUy + wUj; 


12 



with Eqs. (2)-(i|), It follows that 


Px + * vvx + ww^) = 0 


The Euler minlnum principle Is based on elimination of the 
local derivation by using the chain rule 

^ P(j> + ^cj>) = ^ ^22) 

This same equation Is obtained If the equlbllbrlum of forces 
is formed In the y direction or 2 direction. 


6.2 The Pressure Minimum Integral 

That Eq. (22) satisfies the continuity equation can be seen 
If, similar to what we did In section ^.2, we partially Integrat 
over the region of flow and In so doing Eq. (I 3 ) Is taken Into 
account : 


JJJp^ dxdydz 


dxdydz = 0 


Thus the second Integral obviously becomes 
^ Cl)ct)J^ dydz - nf4>c|)^5''^x dxdydz 

’JJJ dxdydz 

Now we can obviously see that 

a<l> 


The second Integral disappears because of Eq. (1), likewise 
the first Integral, since first of all the kinematic flow 
conditions on the boundary of the body has to be satisfied and, 
secondly, the mass flow through the distant boundary must 
disappear. Thus the following simple result holds: 

JJJpcj) dxdydz = 0 


or because of Eq, (22): 

VV(^ + ww(|)) dxdydz = 0 

(23) 


7. Numerical Evaluation /14 

7.1 The Eight-Node Element Cube 

To evaluate Eq. (23) we make use of the finite element 
method. In so doing, the potential is represented in a stepwise 
manner but continuously and without gaps in a three-dimensional 
mesh network. 

A useful approximation for this is the trlllnear element 
cube with the corner point coordinates in the Cartesian auxiliary 
coordinate system C, n, C according to the following table 


i li ^i 


1 

2 

3 

4 

5 

6 


1 


-1 


-1 


1 

1 


-1 


7 -1 

8 1 


-1 -1 


-1 

1 

1 


-1 


-1 


1 

1 


-1 


-1 

-1 

1 

1 

1 

1 
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A function f is then approximated in the cube Including the sides 
by the following equation 


f ■ ^ in (1 ♦Till)' * EiP 

R 

1 


as one can easily convince oneself by trying it out or by 
coefficient determination from the general equation. 


7.2 The General Eight-Node Volume Element /15 

The transition to a volume element with straight edges and 
eight corners is made by successively replacing the function f 
in Eq. (24) by x, y, z and where now the coordinates n, 
and C function as parameters: 

8 

1 

8 

Y =^Yi 
1 

8 

z 

1 

<j) = Gi 

1 

To evaluate Eq. (23) we need the following derivations 

u = <|)x C X 

V = <|)y =(i>g gy T^y +<t>^ £y 

= <t>z 62 Lz 
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In this connection, the local derivations E „ _ and 

x,y,z A,y,z 

If w pose difficulties, since the approximation for x,y and z 
Is not Invert able. The problem can be solved, however. If first 
of all we write down the differentials of all the coordinates. 

dg - 5 y dy + § z dz 

^ y dy + ^ z dz 



dx + 

Ey dy + 

Cz 

dx . Xg 

d^ + 

^T1 * 

"C % 

dy = yg 

dg * 

YT) ^'f\ * 


dz = zg 

dg * 

z^ diq + 



If we now substltue the differentials d-, d^, d_ of the 

t n C 

three last equations by the right sides of the three first 

equations and In each case perform the coefficient comparison 

for dx, dy and dz, we obtain nine equations for the nine unknowns 

tf T* f ,, The time consuming but trivial 

x,y,z x,y,z x,y,z 

solution leads to the following result; 


- y[; z^q ) 

D 

n -C > 

D 

(yg Z7| - y7| zg ) 

D 



(continued on next page) 



Iz 


- y-q ) 

0 

tygxc; - y^xg ) 

D 

txg VTI - XT) yg ) 

0 

with ° - =^1, y^ ) * yg tzi,x^-x-,|2^ ) 

+ z^ (x^i y^ _ y^ X^ ) 


7»3 The Element Pressure Minimum Integral /17 

With the above operations we can now approximately evaluate 
Integral (23) for one element. To do this we use simple point 
Integration by forming the Integrand at point 

g=Tl-C- 0 


and the volume element Is distorted with the Image determinant 
D derived In the previous section: 

dxdydz « D 8 D 


With the exception of an unimportant constant factor, Eq. (23) 
Is then expressed as follows If, for example, we differentiate 
according to the potential for corner 1: 


E 

+ ^^gSz +4>T1 ^z ^4>^£z 


dxdydz «« 





^1-q ^x 

* 

Cx> 

) (Gig gy 



£y> 

^ ^ z 

^iTj^z 

+ Gig 

U>] 

for g.Tj 

a g - 0 
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By rearrangement we obtain a term of the following form for the 
Integral ; 

§(C()| X + (j)7|Y z) 

Finally, If we explicitly differentiate >, we obtain 

J-i 

(X + Y + Z 

1 

7.^ The Global Pressure Minimum Integral /l8 

Now after the pressure minimum integral for an element with 
the pivot point marked ”1” Is evaluated, the construction of the 
entire volume integral no longer presents any difficulties. 

All that Is necessary Is for the contributions of all of the 
elements which contain the corner ”1” to be counted together. 

So In the free field the contributions of eight elements, which 
have one and the same corner in common, are to be counted 
together, as one can easily Imagine. Boundaries present no 
problems. Here, of course. In most cases less than eight elements 
are Involved In the integrand. In general, four elements are 
Involved. 


7.5 Relaxation Method 

So far we have treated Eq. (23) for the case of potential 
derivation at a selected mesh corner point. In order to solve 
Eq. (23) in a global association of elements, we start with the 
idea that the pressure minimum Integral must disappear for 
each Intersection of the mesh. Thus, from Eq. (23) we obtain Just 
as many equations as unknown potential given values. 

Prom this we can derive in a classical manner a relaxation 
method; 


18 


Up to a certain distance the flow field Is completely 
covered with volume elements. The distant field solution of 
the potential as per Eqs. (16), (20) and (21) Is specified on 
the distant field border. The trldlagonal matrix for the 
unknown potential values is set up along reasonably selected 
lines . 


All contributions not stemming from this line are put on the 
right side. The resulting system of equations is solved directly. 

If we now move from line to line, we obtain an iterative 
algorithm with which the rewest (> and p values are always used. /19 
In so doing, the flow field Is repeatedly traversed until the 
potential values no longer change significantly. 

Caution Is necessary in determining the potential values 
on the top suid bottom of the vortex layer. In this computer 
program, numerical stability was achieved as follows; 


We define a mean potential as follows; 


<!> 


m 




and determine the potential values on the top and bottom of the 

trailing edge of the wing ^uHK’ constant 

potential Jump ■ ^oHK “ ^uHK sections y ■ constant are 

assumed to be known. Now In order to calculate the potential 

values on both sides of the vortex layer, we first of all solve 

for 6 and then form 
m 


4^0 

4>u 


A 

9m 

1 

9tn “ 2 
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8« Results 


NACA~0012 Wing 

Fop the first test sample we used as a basis for the 
calculation a simple swept back wing with linear warping and 
a constant profile. 

Pig. 1 shows the outer boundary belonging to each y section 
on which the distant field potential according to Eqs. (16), 

(20) and (21) Is specified. The Influences of the other half of 
the wing on the computation fields shown are also taken Into 
consideration. 

The wing Is centered on an impenetrable wall so that even 
flows with small slldesllp angles can be calculated with a very 
simplified fuselage effect. Beyond the wlngtlp the calculating 
network Is expanded by another three sections which contain no 
solid contours. The vortex layer leaves the trailing edge in 
the direction of the local bisecting lines. 

Figs. 2-l6 show the pressure distributions. The points 
above the profile indicate the elements in which supersonic 
velocity prevails. The mark on the C axis represents the 

tr 

critical pressure coefficient for the total oncoming flow Mach 
number. 

Pig. l8 requires more detailed explanation. The top side 
pressure distribution acts somewhat unevenly, since the shock 
front ;)umps along the wingspan a few times between the mesh 
coordinates . 

In this as In all other shock capturing methods, the shock 
front can coincide only with network coordinates so that , In 
general, the Ranklne-Hugonlot shock equation and the angle 
correlations in front of and behind the shock cannot be satisfied. 


If one wanted to Improve this, the method would be very com- 
plicated and would use an enormous amount of computer time. We /21 
can manage with this method by approximately anticipating all 
possible shock formations by means of a mesh construction per- 
pendicular to the body. Even this is very difficult. It is 
simpler to concentrate the network in such a way that shocks can 
be geometrically approximated in steps. Because of a lack of 
computer time, no studies on this could be carried out. 

The problem is most easily coped with by intensifying the 
artificial viscosity via the specifiable viscosity parameters. 

Then the result becomes almost Independent of the ’’falde” local 
network geometry in the vicinity of the shock. 

Figs. 20-22 are auxiliary diagrams showing the coefficient 
distributions with respect to the local chord length. In these 
diagrams the pitching moment is formed around the forward most 
point of the wing. Using this, together with the lift distribution 
and position of the pressure point, we can convert the moment to 
any point of impact. 

The iterative convergence is considerably accelerated through 
successive network divisions with subsequent potential inter- 
polation. So for this example the result was obtained after eight 
minutes of computer time on the cnetral unit of an IBM 370/168. 

This required only 12 iterations per network (three altogether). 

The computer time could be reduced even more considerably if 
greater access storage were available. In the absence of storage 
capacity, the network geometry and the transformation matrix for 
each element must now be calculated anew for each iteration. 

However, this is still cheaper than reading off this once deter- 
mined data from an external and thus slow storage. 
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This is a design study for the wing of a modern transport 
plane. Again the same data are supplied as In the previous 
example . 








In particular. Fig, 36 shows that obviously a well balanced 
upper wing surface has been achieved, for the shock Intensities 
on the outer wing are very small. At about two-thirds of the 
chord length (y “ 11.75, Fig. 31) the pressure distribution of 
the basic profile Is satisfactorily represented. 


This example with 12 wing sections for each 32 profile 
points required the solution of ^i ,320 unknown potential values. 
The computer time was five minutes on the central computer of the 
IBM 370/168. 


9 . Summary /23 

This paper demonstrates the operability of a new finite 
volume method for calculating shock-affected flow around a wing. 
Because of the chronic lack of computer time no calculations could 
be carried out on the dense mesh network provided by the program 
which would certainly have produced more accurate results. 
Nevertheless, on the basic of earlier experiments. It is worth- 
while to extend the program to wing-fuselage combinations in w 
which, because of the complicated boundary geometry, the 
advantages of the finite volume technique actually first come to 
bear* For this taks, the usual difference methods are suitable 
ortly to a limited extent, since they all require a strictly 
perpendicular network which leads to a considerable amount of 
Interpolation on solid surfaces. For this reason and because of 
the quicker Iterative convergence, methods such as this one 
should take on considerable Importance in the future. 
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